#pragma once
#include <vector>
#include <iostream>
#include "MVoperator.h"

using namespace std;

MatrixVectorOperator MV;

class Function {
public:
    virtual vector<double> operator()(vector<double> x) = 0;
};

class InterpolationOperator {
public:
    virtual vector<double> interpolate(vector<double> &coarse) = 0;
};

class LinearInterpolationOperator : public InterpolationOperator {
public:
    vector<double> interpolate(vector<double> &coarse) {
        int n = coarse.size() + 1;
        vector<double> fine(2 * n - 1);
	fine = MV.M_times_V(MV.get_Iu(n),coarse);

        return fine;
    }
};

class QuadraticInterpolationOperator : public InterpolationOperator {
public:
    // 将粗网格coarse插值到长度为2n-1的细网格fine上
    vector<double> interpolate(vector<double> &coarse) {
	int n = coarse.size()+1;
	int n_h = 2 * n - 1;
	vector<double> fine(n_h);
	double h = 1.0 / (n-1);
        double H = 2.0 * h;
        fine[0] = 0.25 * coarse[n-2] + 0.75 * coarse[1];

        // 对相邻节点之间的值进行二次插值
        for (int i = 2; i <= n_h-2; i += 2) {
            double a = (coarse[i/2-1] - 2*coarse[i/2] + coarse[i/2+1]) / (h*h);
            double b = (coarse[i/2+1] - coarse[i/2]) / h - 0.5 * h * (2*a + coarse[i/2+1] - coarse[i/2]);
            double c = coarse[i/2] - 0.25 * h * h * a - 0.5 * h * b;
            double x_i = (i-1) / 2.0 * H;
            fine[i] = a*(x_i-H)*(x_i-H) + b*(x_i-H) + c; 
            fine[i+1] = a*x_i*x_i + b*x_i + c; 
        }

        fine[n_h-1] = 0.75 * coarse[n-2] + 0.25 * coarse[n-3];
	return fine;
    };
};

class RestrictionOperator {
public:
    virtual vector<double> restrict(vector<double> &fine) = 0;
};

class FullWeightingRestrictionOperator : public RestrictionOperator {
public:
    vector<double> restrict(vector<double> &fine) {
        int n = (fine.size() + 1) / 2;
	vector<double> coarse(n-1);
	coarse = MV.M_times_V(MV.get_Id(2*n),fine);
        return coarse;
    }
};

class InjectionRestrictionOperator : public RestrictionOperator {
public:
    vector<double> restrict(vector<double> &fine) {
        int n = (fine.size() + 1) / 2;
        vector<double> coarse(n - 1);

        coarse[0] = fine[0];
        coarse[n - 2] = fine[2 * n - 2];
        for (int i=1; i < n - 2; i++) {
            coarse[i] = fine[2 * i];
        }

        return coarse;
    }
};
